In [1]:
import numpy
from matplotlib import pyplot
%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
from JSAnimation.IPython_display import display_animation
from matplotlib import animation
In [2]:
# Constants
nx = 81
dx = 0.25
dt = 0.0002
gamma = 1.4 # adiabatic index
# initial data
RHO_L = 1. #kg/m^3
U_L = 0. # m/s
P_L = 100000. # N/m^2
# right initial data
RHO_R = 0.125 # kg/m^3
U_R = 0. # m/s
P_R = 10000. # N/m^2
XMIN=-10 # m
XMAX = 10 # m
X_0 = 0 # m, point where the membrane is at t = 0
# The left primitive vector
V_PL = numpy.array([RHO_L,U_L,P_L])
# The right primitive vector
V_PR = numpy.array([RHO_R,U_R,P_R])
We define two vectors, a primitive vector v_primitive and a conserved vector v_conserved. The primitive vector contains rho, u, and p, while the conserved vector contains rho,rho u,rho e_T, where e_T is the total energy given by the equation of state.
First we need to define maps from the pimitive to the conserved variables, and vice versa.
In [3]:
def get_p(v_conserved):
"""
Returns the pressure as a function of the
conserved vector.
"""
v = v_conserved # for convenience
return (gamma-1)*(v[2]-0.5*(v[1]**2)/v[0])
In [4]:
def get_e_T(v_primitive):
"""
Returns the energy as a function of the primitive vector
"""
rho = v_primitive[0]
u = v_primitive[1]
p = v_primitive[2]
e = p / ((gamma - 1)*rho)
e_T = e + 0.5 * (u**2)
return e_T
In [5]:
def get_conserved(v_primitive):
"""
Returns the conserved vector
based on the primitive vector
"""
rho = v_primitive[0]
u = v_primitive[1]
p = v_primitive[2]
e_T = get_e_T(v_primitive)
return numpy.array([rho,rho*u,rho*e_T])
In [6]:
def get_primitive(v_conserved):
"""
Returns the vector of primitive variables as a
function of the conserved variables.
"""
rho = v_conserved[0]
u = v_conserved[1]/rho
e_T = v_conserved[2]/rho
p = get_p(v_conserved)
return numpy.array([rho,u,p])
We also compute the flux
In [7]:
def flux_of_primitive(v_primitive):
"""
Computes the flux as a function of the primitive vector.
"""
rho = v_primitive[0]
u = v_primitive[1]
p = v_primitive[2]
e_T = get_e_T(v_primitive)
return numpy.array([rho*u,
rho*u*u+p,
u*(rho*e_T+p)])
In [8]:
def compute_flux(v_conserved):
"""
Computes the flux as a function of the
conserved variables.
"""
v = v_conserved
out = numpy.empty_like(v)
out[0] = v[1]
out[1] = ((v[1]**2)/v[0]
+ (gamma - 1)*(v[2] - 0.5*v[1]*v[1]/v[0]))
out[2] = (v[2]
+ (gamma-1)*(v[2] - 0.5*v[1]*v[1]/v[0]))*(v[1]/v[0])
return out
In [9]:
x = numpy.linspace(XMIN,XMAX,nx)
In [10]:
def get_v0_primitive(x):
"""
Returns the initial data as a function of x in
primitive form.
"""
if x < X_0:
return V_PL
else:
return V_PR
In [11]:
def richtmyer_step(v0,dt,dx):
"""
steps v0 forward one step in time
and returns the new solution
"""
v_out = numpy.empty_like(v0)
v_half = numpy.empty_like(v0)
flux = numpy.array(map(compute_flux,v0))
v_half[:-1] = 0.5 * ((v0[:-1] + v0[1:])
- (dt/dx)*(flux[1:]-flux[:-1]))
v_half[-1] = v0[-1]
flux_half = numpy.array(map(compute_flux,v_half))
v_out[1:] = (v0[1:]
- (dt/dx)*(flux_half[1:]-flux_half[:-1]))
v_out[0] = v0[0] # dirichlet boundary conditions
v_out[-1] = v0[-1]
return v_out
In [12]:
def richtmyer(v0,nt,dt,dx):
"""
Computes the solution for the conserved variables
with a Richtmyer scheme.
INPUTS
-------
v0 is the initial conserved vector
nt is the number of times
dt is the change in time per time step
dx is the change in space per grid point
OUTPUTS
--------
times, a one-dimensional array of hte times
(initial time is assumed to be zero)
v_n, a matrix where each row is the v_conserved
solution at a given time.
"""
# initialize our arrays
v_n = numpy.zeros((nt,
v0.shape[0],
v0.shape[1]))
times = numpy.zeros(nt)
v_n[0] = v0.copy()
times[0] = 0.
# main loop
for i in range(1,nt):
v_n[i] = richtmyer_step(v_n[i-1],dt,dx)
times[i] = times[i-1] + dt
# output
return times, v_n
In [17]:
v0_primitive = numpy.array(map(get_v0_primitive,x))
v0_conserved = numpy.array(map(get_conserved,v0_primitive))
nt = int(.01/dt)+1
nt
Out[17]:
In [18]:
times,v_n = richtmyer(v0_conserved,nt,dt,dx)
In [19]:
v_np = numpy.empty_like(v_n)
flux_np = numpy.empty_like(v_n)
for i in range(nt):
v_np[i] = numpy.array(map(get_primitive,v_n[i]))
flux_np[i] = numpy.array(map(compute_flux,v_n[i]))
In [20]:
rho_n = v_np[...,...,0]
u_n = v_np[...,...,1]
p_n = v_np[...,...,2]
In [23]:
pyplot.plot(x,p_n[50])
Out[23]:
In [26]:
numpy.where(x==2.5)
Out[26]:
In [25]:
x[50],times[-1]
Out[25]:
In [30]:
v_np[-1,50]
Out[30]:
In [ ]: